SUBROUTINE G_MOT_R
!-----------------------------------------------------------------------
!
! Definition von G:  Bewegungs- bzw. Impulsgleichung
!

      use primvar,  only : G, X, XZ, XA, MD, MUr, MUp, MR, MH, ME, Gnorm
      use physco,   only : sqrt2pi, pi, grav, offset, clight, z1, z2, z4, z12
      use config,   only : np, BCflag
      use global,   only : tst, relax
      use geomvar,  only : V_flux, V_vol, V_volA, V_volZ, S_volZ
      use matvar,   only : Pgas0Z, OPArosZ
      use viscvar,  only : Uterm_r, muQ_r
      use advecvar, only : momr_adv
      use zvar,     only : H_pZ
      use RBvar,    only : M_star


      implicit none

      integer       :: i
      real(kind=8) :: inflow_ur
      logical, save :: init = .true.

      if (init) then
         if (.not. relax) then
            write(66,"(a)") "G_mot_r.f90:   Innenrand: Gl. loesen ;  Aussenrand: Gl. loesen"
         else
            write(66,"(a)") "G_mot_r.f90:   Zum Relaxieren: Innenrand: u_r = 0 ;  Aussenrand: u_r = 0"
         end if
         init = .false.
      end if


!-----------------------------------------------------------------------
!    Randbedingungen
!-----------------------------------------------------------------------

      inflow_ur = 4.e+03

   ! innere Pseudozellen: i=1 & i=2
      G(MUr,1)        = X(MUr,1)
      Gnorm(MUr,1)    = max( offset, abs( X(MUr,1) ) )
      G(MUr,2)        = X(MUr,2)
      Gnorm(MUr,2)    = max( offset, abs( X(MUr,2) ) )

   ! aeussere Pseudozellen: i=np, i=np-1
      G(MUr,np)       = X(MUr,np) + inflow_ur 
      Gnorm(MUr,np)   = max( offset, abs( X(MUr,np) + inflow_ur  ) )
      G(MUr,np-1)     = X(MUr,np-1) + inflow_ur
      Gnorm(MUr,np-1) = max( offset, abs( X(MUr,np-1) + inflow_ur ) )
      G(MUr,np-2)     = X(MUr,np-2) + inflow_ur
      Gnorm(MUr,np-2) = max( offset, abs( X(MUr,np-2) + inflow_ur ) )


!-----------------------------------------------------------------------
!    Restlicher Bereich
!-----------------------------------------------------------------------

! VORSICHT: X(MD,np-2) ist schon die Dichte-Randbedingung!
      do i=3,np-3

         G(MUr,i) = V_vol(i) * X(MUr,i) * X(MD,i) - V_volA(i) * XA(MUr,i) * XA(MD,i)                                 & ! temporal diff.
                  + V_flux(i) * momr_adv(i) - V_flux(i-1) * momr_adv(i-1)                                            & ! advection
                  - z12 * ( XZ(MUp,i)**2 * S_volZ(i) + XZ(MUp,i-1)**2 * S_volZ(i-1) ) / XZ(MR,i) * XZ(MD,i) * tst    & ! zentrifugal force
                  + z2*pi * XZ(MR,i) * sqrt2pi * ( H_pZ(i) * Pgas0Z(i) - H_pZ(i-1) * Pgas0Z(i-1) ) * tst             & ! gas pressure force
                  + grav * M_star / XZ(MR,i)**2 * XZ(MD,i) * V_volZ(i) * tst                                         & ! grav. accel.
                  - z2*pi / clight * XZ(MH,i) *                                                                      & ! rad. pressure force
                          ( OPArosZ(i)*XZ(MD,i)*S_volZ(i) + OPArosZ(i-1)*XZ(MD,i-1)*S_volZ(i-1) ) * tst              & ! ...
                  - z12*pi / XZ(MR,i) *  (   muQ_r(i)   * ( XZ(MR,i)**2 + XZ(MR,i+1)**2 ) * Uterm_r(i)   -           & ! viscosity
                                             muQ_r(i-1) * ( XZ(MR,i)**2 + XZ(MR,i-1)**2 ) * Uterm_r(i-1)   ) * tst     ! ...


         Gnorm(MUr,i) = max(                                                                                         &
               abs( V_vol(i) * X(MUr,i) * X(MD,i) ),                                                                 & ! temporal diff.
               abs( V_flux(i) * momr_adv(i) ), abs( V_flux(i-1) * momr_adv(i-1) ),                                   & ! advection
               abs( z12 * ( XZ(MUp,i)**2 * S_volZ(i) + XZ(MUp,i-1)**2 * S_volZ(i-1) ) / XZ(MR,i) * XZ(MD,i) * tst ), & ! zentrifugal force
               abs( z2*pi * XZ(MR,i) * sqrt2pi * ( H_pZ(i) * Pgas0Z(i) - H_pZ(i-1) * Pgas0Z(i-1) ) * tst ),          & ! gas pressure force
               abs( grav * M_star / XZ(MR,i)**2 * XZ(MD,i) * V_volZ(i) * tst ),                                      & ! grav. accel.
               abs( z2*pi / clight * XZ(MH,i) *                                                                      & ! rad. pressure force
                          ( OPArosZ(i)*XZ(MD,i)*S_volZ(i) + OPArosZ(i-1)*XZ(MD,i-1)*S_volZ(i-1) ) * tst ),           & ! ...
               abs( z12*pi / XZ(MR,i) *  (   muQ_r(i)   * ( XZ(MR,i)**2 + XZ(MR,i+1)**2 ) * Uterm_r(i)   -           & ! viscosity
                                             muQ_r(i-1) * ( XZ(MR,i)**2 + XZ(MR,i-1)**2 ) * Uterm_r(i-1)   ) * tst ) & ! ...
                          )

      end do


!-----------------------------------------------------------------------
!    Geschlossene Raender fuers Relaxieren
!-----------------------------------------------------------------------

! Ueberschreibt einfach Randwerte aus obiger Schleife
   ! aeussere Randbedingung: i=np-2                                                                   !debugging
      G(MUr,np-2)     = X(MUr,np-2)  + inflow_ur            !    fixer Randwert - u_r = 0             !debugging
    
      Gnorm(MUr,np-2) = MAX( offset, ABS(X(MUr,np-2) + inflow_ur ) )                                !debugging
!   WRITE(*,*) 'Ur am Außenrand = ', X(MUr,np-2), 'gesetzter Wert = ' ,inflow_ur
    
      if (relax) then
      ! innere Randbedingung: i=3
         G(MUr,3)     = X(MUr,3)                              !    fixer Randwert - u_r = 0
         Gnorm(MUr,3) = MAX( offset, ABS(X(MUr,3) ) )
      ! aeussere Randbedingung: i=np-2
         G(MUr,np-2)     = X(MUr,np-2)                        !    fixer Randwert - u_r = 0
         Gnorm(MUr,np-2) = MAX( offset, ABS(X(MUr,np-2) ) )
      end if


END SUBROUTINE G_MOT_R

